

# This was run using a Macbook Pro with an Apple M1 Pro chip, using R version 4.2.2 in RStudio 2025.05.1+513.
# Requires the installation of tjbal (https://yiqingxu.org/packages/tjbal/tjbal.html).


rm(list=ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

library(tjbal)
library(haven)


          #################
          ##### TJBAL #####
          #################

            ############
            # Case T=G #
            ############

#Load data
DiDdata <- read_dta("../dta/did_rep.dta")

# Subset the data
DiDdata<-DiDdata[!(DiDdata$sch_A == 1 | DiDdata$sch_B == 1),]
DiDdata<-DiDdata[is.na(DiDdata$drummond),]
DiDdata<-DiDdata[!(DiDdata$year==1830 | DiDdata$year==1831 | DiDdata$year==1832),]

# Define variables
DiDdata$post<-0
DiDdata$post[DiDdata$year>1832] <-1
DiDdata$treatgroup<-0
DiDdata$treatgroup[DiDdata$new == 1] <- 1

# Define treatment
DiDdata$treated <- DiDdata$post * DiDdata$treatgroup


# Mean balancing
out.did <- tjbal(ln_petitions ~ treated, data=DiDdata, 
                 index=c("ID","year"), demean=FALSE, estimator="mean", vce="jackknife")
print(out.did)

# Plots
jpeg(filename="figure4a.jpeg")
plot(out.did, type="counterfactual", count=FALSE)
dev.off()

jpeg(filename="figure4b.jpeg")
plot(out.did, ylab="ATT", ylim=c(-1.5,1))
dev.off()


          ############
          # Case T=L #
          ############

# Load data
DiDdata <- read_dta("../dta/did_rep.dta")

# Subset the data
DiDdata<-DiDdata[!(DiDdata$new==1),]
DiDdata<-DiDdata[DiDdata$sch_A==1 | DiDdata$sch_B==1 | !is.na(DiDdata$drummond),]
DiDdata<-DiDdata[!(DiDdata$year==1830 | DiDdata$year==1831 | DiDdata$year==1832),]

# Define variables
DiDdata$post<-0
DiDdata$post[DiDdata$year>1832] <-1
DiDdata$treatgroup<-0
DiDdata$treatgroup[DiDdata$sch_A == 1] <- 1

# Define treatment
DiDdata$treated <- DiDdata$post * DiDdata$treatgroup

# Mean balancing
out.did <- tjbal(ln_petitions ~ treated, data=DiDdata, index=c("ID","year"), demean = FALSE, estimator="mean", vce="jackknife")
print(out.did)

# Plots
jpeg(filename="figure4c.jpeg")
plot(out.did, type="counterfactual", count=FALSE)
dev.off()

jpeg(filename="figure4d.jpeg")
plot(out.did, ylab="ATT", ylim=c(-1,0.5))
dev.off()














